#pragma once

#include <complex>
#include <vector>
#include "nnec/nnec_errors.h"

using namespace std;

namespace NNEC
{
    template<typename T>
    struct is_complex_or_floating_point : is_floating_point<T>{};

    template<typename T>
    struct is_complex_or_floating_point<complex<T>> : is_floating_point<T>{};

    class INTEGRAND {
    public:
        double (*function) (double x, void *para){nullptr};
        void * parameter{nullptr};
    };

    int integration_romberg(INTEGRAND * integrand, double a, double b, int n,
                            double epsabs, double epsrel, double * result, int *eval_count);

    vector<double> vector_add(const vector<double> & a, const vector<double> &b);
    vector<double> vector_sub(const vector<double> & a, const vector<double> &b);
    vector<double> vector_scale_mul(double a, const vector<double> & v);
    double vector_dot_mul(const vector<double> & a, const vector<double> & b);
    vector<double> vector_cross_mul(const vector<double> & a, const vector<double> & b);
    double vector_norm(const vector<double> & a);
    vector<double> vector_normalize(const vector<double> & a);

    template<typename T>
    int lu_decomposition(T *A, int n, int *permutation, int *d)
    {
        static_assert(is_complex_or_floating_point<T>::value, "The template type is not supported.");

        if(n < 1 || A == nullptr || permutation == nullptr || d == nullptr) {
            return NNEC_ERROR_INVALID;
        }

        unique_ptr<double[]> vmax(new double[n]);
        if(!vmax) {
            return NNEC_ERROR_ALLOCATION;
        }

        int i_max;
        *d = 1;

        for(int i = 0; i < n; i++) {
            double row_max = 0.0;
            for(int j = 0; j < n; j++) {
                double tmp = abs(A[i * n + j]);
                if(tmp > row_max) {
                    row_max = tmp;
                }
            }
            if(row_max == 0.0) {
                return NNEC_ERROR_SINGULARITY;
            }
            vmax[i] = 1.0 / row_max;
        }

        for(int j = 0; j < n; j++) {
            for(int i = 0; i < j; i++) {
                T sum = A[i * n + j];
                for(int k = 0; k < i; k++) {
                    sum = sum - A[i * n + k] * A[k * n + j];
                }
                A[i * n + j] = sum;
            }

            double row_max = 0.0;
            for(int i = j; i < n; i++) {
                T sum = A[i * n + j];
                for(int k = 0; k < j; k++) {
                    sum = sum - A[i * n + k] * A[k * n + j];
                }
                A[i * n + j] = sum;
                double tmp = vmax[i] * abs(sum);
                if(tmp >= row_max) {
                    row_max = tmp;
                    i_max = i;
                }
            }

            if(j != i_max) {
                for(int k = 0; k < n; k++) {
                    T tmp = A[i_max * n + k];
                    A[i_max * n + k] = A[j * n + k];
                    A[j * n + k] = tmp;
                }
                *d = -(*d);
                vmax[i_max] = vmax[j];
            }
            permutation[j] = i_max;
            if(abs(A[j * n + j]) == 0.0) {
                return NNEC_ERROR_SINGULARITY;
            }
            if(j != n - 1) {
                T tmp = 1.0 / A[j * n + j];
                for(int i = j + 1; i < n; i++) {
                    A[i * n + j] = A[i * n + j] * tmp;
                }
            }
        }
        return NNEC_SUCCESS;
    }

    template<typename T>
    int lu_backward_substitution(T *A, int n, int const *permutation, T *b)
    {
        static_assert(is_complex_or_floating_point<T>::value, "The template type is not supported.");

        if(n < 1 || A == nullptr || permutation == nullptr || b == nullptr) {
            return NNEC_ERROR_INVALID;
        }

        int nv = -1;
        for(int i = 0; i < n; i++) {
            int ip = permutation[i];
            T sum = b[ip];
            b[ip] = b[i];
            if(nv >= 0) {
                for(int j = nv; j <= i - 1; j++) {
                    sum = sum - A[i * n + j] * b[j];
                }
            }
            else if(abs(sum) > 0.0) {
                nv = i;
            }
            b[i] = sum;
        }

        for(int i = n - 1; i >= 0; i--) {
            T sum = b[i];
            for(int j = i + 1; j < n; j++) {
                sum = sum - A[i * n + j] * b[j];
            }
            b[i] = sum / A[i * n + i];
        }

        return NNEC_SUCCESS;
    }
}
